home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / rzextr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  859 b   |  51 lines
  1. /*
  2. ### rational extrapolation ###
  3. */
  4.  
  5. extern double **ddd, *xxx;
  6.  
  7. void rzextr(iest,xest,yest,yz,dy,nv,nuse)
  8. int iest,nv,nuse;
  9. double xest,yest[],yz[],dy[];
  10. {
  11.     int  m1,k,j;
  12.     double yy,v,ddy,c,b1,b,*fx,*dvector();
  13.  
  14.     fx = dvector(0,nuse-1);
  15.     xxx[iest] = xest;
  16.     if(iest == 0)
  17.         for(j=0;j<nv;j++) {
  18.             yz[j] = yest[j];
  19.             ddd[j][0] = yest[j];
  20.             dy[j] = yest[j];
  21.         }
  22.     else {
  23.         m1 = (iest < nuse-1 ? iest : nuse-1 );
  24.         for(k=0;k<m1;k++)
  25.             fx[k+1] = xxx[iest-k-1] / xest;
  26.         for(j=0;j<nv;j++){
  27.             yy = yest[j];
  28.             v = ddd[j][0];
  29.             c = yy;
  30.             ddd[j][0] = yy;
  31.             for(k=1;k<m1+1;k++){
  32.                 b1 = fx[k] * v;
  33.                 b = b1 - c;
  34.                 if (b) {
  35.                     b = (c - v) / b;
  36.                     ddy = c * b;
  37.                     c = b1 * b;
  38.                 }
  39.                 else 
  40.                     ddy = v;
  41.                 if(k != m1) v = ddd[j][k];
  42.                 ddd[j][k] = ddy;
  43.                 yy += ddy;
  44.             }
  45.             dy[j] = ddy;
  46.             yz[j] = yy;
  47.         }
  48.     }
  49.     (void) free_dvector(fx,0,nuse-1);
  50. }
  51.